
d=load("6.28.all")
place=d(:,1)
lati0=d(:,2)*1e-5+30.74
longi0=d(:,3)*1e-5+103.92
latitude=d(:,4)*1e-5+30.74
longitude=d(:,5)*1e-5+103.92

radiusEquatorial = 6378.1370*1000
radiusPolar = 6356.7523*1000

function rad=toRadians(degrees)
rad=degrees*pi/180
endfunction

function root=radius(x,y)
root=sqrt(power(x, 2) + power(y, 2))
endfunction

function r=earthRadius(latitude)
theta = toRadians(latitude)
a=6378.1370*1000 
b=6356.7523*1000
m = radius(power(a, 2) * cos(theta), power(b, 2) * sin(theta))
n= radius(a * cos(theta), b * sin(theta))
r=m/n
endfunction

function dis=distance(lat1,long1,lat2,long2,avg)
dis=2* earthRadius(lat1) * asin(sqrt(power(sin((toRadians(lat2-lat1))/2),2)+cos(toRadians(lat1))*cos(toRadians(lat2))*power(sin((toRadians(long2-long1))/2),2)))*cos(theta(lat1,long1,lat2,long2)-avg)
endfunction

function th=theta(lat1,long1,lat2,long2)
    th=atan2(toRadians(lat2 - lat1), toRadians(long2 - long1) * cos(lat1))
endfunction

#與初始節點相差
#args=[latitude longitude zeros(length(latitude),1)+latitude(1) zeros(length(longitude),1)+longitude(1)]
args=[lati0	longi0 latitude longitude]

thetaAll=arrayfun(@theta,args(:,1),args(:,2),args(:,3),args(:,4))
avg=mean(thetaAll)
result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4),zeros(length(args),1)+avg)

#相鄰兩點相差
#lati=latitude(2:37)-latitude(1:36)
#longti=longitude(2:37)-longitude(1:36)

#args=[lati longti zeros(36,1) zeros(36,1)]
#result=arrayfun(@distance,args(:,1),args(:,2),args(:,3),args(:,4))
#result=abs(result)

plot(place,result,place,place*2.5-1.25)
result=result-place*2.5+1.25
x=0:0.1:10
result=abs(result)
m=mean(result)
s=std(result)
d=power(s,2)
f=(erf((x-m)/(s*sqrt(2)))+1)/2
p=exp(power(x-m,2)/(-2*d))/(s*sqrt(2*pi))
figure
plot(x,f)
#figure
#plot(x,p)
